/* Author : 
 *    Daohang Shi 
 *    dshi7@wisc.edu 
 *    Fall  2013  
 *    */

#include "Parser_cplex.h" 
//---------------------------------------------------------------------------
int Parser_cplex::getCplex ()
{	
  int status = 0;
  srand (time (NULL));
  while (1) 
  { 	                         
    lpEnv = CPXopenCPLEX (&status); 
    if (lpEnv != NULL)
      break;
    else
    {
      char errmsg[1024];
      fprintf (stderr, "Could not open CPLEX environment.\n");
      CPXgeterrorstring (lpEnv, status, errmsg);
      fprintf (stderr, "%s", errmsg);
      goto TERMINATE;
      int time = 450 + rand()%300;
      fprintf (stderr, "Wait %d seconds!!\n", time);
      Sleep (time);
    }
  }	      
  srand (0);
#ifndef CPLEX65	
  //** Turn on output to the screen **
  status = CPXsetintparam (lpEnv, CPX_PARAM_SCRIND, CPX_OFF);
  if (status) 
  {
    fprintf (stderr, 
        "Failure to turn on screen indicator, error %d.\n", status);
    goto TERMINATE;
  }	                
  //** Turn on data checking **
  status = CPXsetintparam (lpEnv, CPX_PARAM_DATACHECK, CPX_OFF);
  if (status) 
  {
    fprintf (stderr, 
        "Failure to turn on data checking, error %d.\n", status);
    goto TERMINATE;
  }
#endif	
  //* Create the problem **
  lp = CPXcreateprob (lpEnv, &status, "GR");	
  if (lp == NULL) 
  {
    fprintf (stderr, "Failed to create LP.\n");
    goto TERMINATE;
  }	     

TERMINATE:	
  return (status);
}

//---------------------------------------------------------------------------
int Parser_cplex::CGR_CreateProblemLP ()
{	

  cout<<endl<<"create LP problem"<<endl << endl;
  int status = 0;
  CPXchgobjsen ( lpEnv, lp, CPX_MIN );

  for ( size_t i=0; i < cgr.bnets.size(); i++ ) 
    cgr.x_var_size += cgr.bnets[i].candidate_routes.size();

  cout << "Total number of candidate net routes = " << cgr.x_var_size << endl;

  status = CGR_AddNewRows (lpEnv, lp);
  assert ( !status );

  status = CGR_AddGridEdgeColumns (lpEnv, lp);
  assert ( !status );

  status = CGR_AddCoarseNetColumns (lpEnv, lp);
  assert ( !status );

  status = CGR_SolveLP (lpEnv, lp);
  assert ( !status );

TERMINATE:	

  return (status);
}
//---------------------------------------------------------------------------
int Parser_cplex::DeleteProblemLP ()
{         
  int status;

  //** Free up the problem as allocated neighbors[3] CPXcreateprob, if necessary **
  if (lp != NULL ) 
  {
    status = CPXfreeprob (lpEnv, &lp);
    if (status) 
      fprintf (stderr, "CPXfreeprob failed, error code %d.\n", status);
  }

  //** Free up the CPLEX environment, if necessary **
  if (lpEnv != NULL) 
  {
    status = CPXcloseCPLEX (&lpEnv);

    if (status > 0) 
    {
      char  errmsg[1024];
      fprintf (stderr, "Could not close CPLEX environment.\n");
      CPXgeterrorstring (lpEnv, status, errmsg);
      fprintf (stderr, "%s", errmsg);
    }
  }	

  //delete [] dualSol;
  dualSol = NULL;
  Free_and_Null ((char **) &cmatbeg);
  Free_and_Null ((char **) &cmatind);
  Free_and_Null ((char **) &cmatval);
  Free_and_Null ((char **) &lb);
  Free_and_Null ((char **) &ub);
  Free_and_Null ((char **) &primalSol);
  Free_and_Null ((char **) &reducedCost);
  Free_and_Null ((char **) &obj);
  Free_and_Null ((char **) &dualSol);

  lpEnv = NULL;
  lp = NULL;
  primalSol = NULL;
  dualSol = NULL;    
  reducedCost = NULL;
  ilpSol = NULL;

  return  status;
}

/* 
 * This function is the extremely important one that solves ILP and also uses the solution back to routing problem.
 * Besides the solution obtained by CPLEX, it is more important to clarify how we use the primal solution and dual solution.
 *  1.  PrimSol Obj is the objective value we are always trying to reduce. And it is reduced or unchanged after each iteration.
 *  2.  The primal solution is a vector contains te and xit in order. The first part of te add up to objective value, and the second part of xit is the decision variables. For each net, pick the entry with largest value as the base route for the following iteration. 
 *  3.  The dual solution are used to define the edge weights for next iteration. 
 *  */
int Parser_cplex::CGR_SolveLP (CPXENVptr env, CPXLPptr lp)
{
  int status;

  cout<<endl<<"solving LP problem"<<endl;

  /* compute the total number of candidate routes */
  cgr.x_var_size = 0;
  for ( vector<CoarseNet>::iterator itr_net = cgr.bnets.begin(); itr_net != cgr.bnets.end(); itr_net ++ ) {
    cgr.x_var_size += (*itr_net).candidate_routes.size();
  }

  status = CPXsetdblparam (lpEnv, CPX_PARAM_TILIM, 10);
  if (status) {
    cerr << "Cannot set time limitation : " << status << endl;
    return  -1;
  }			

  status = CPXsetintparam (lpEnv, CPX_PARAM_LPMETHOD , CPX_ALG_BARRIER);
  if (status) {
    cerr << "Cannot set optimization method : " << status << endl;
    return  -1;
  }

  status = CPXlpopt (env, lp);
  if (status) {
    cerr << "Cannot optimize LP problem : " << status << endl;
    return  -1;
  }

  primalSol = new double [cgr.E_id.size() + cgr.x_var_size];

  dualSol = new double [cgr.bnets.size() + cgr.E_id.size() * cgr.jsize];

  numcols = CPXgetnumcols ( lpEnv, lp );
  numrows = CPXgetnumrows ( lpEnv, lp );

  cout << "In this iteratoin : " << endl;
  cout << "\tnumber of congested edges = " << cgr.congested_edges.size() << endl;
  cout << "\tnumcols = " << numcols << endl;
  cout << "\tnumrows = " << numrows << endl;
  cout << "\tobjValue = " << objValue << endl;
  status = CPXsolution (env, lp, &solStat, &objValue, primalSol, dualSol, NULL, NULL );

  bool  debug_mode = false;
  double obj;

  /* update net route decision variables */
  int primal_index = cgr.E_id.size();

  for ( size_t i=0; i< cgr.bnets.size(); i++ ) {
    //  pick the route with largest decision variable
    int index_of_max_decision_variable = -1;
    if ( debug_mode==true )
      cout << "bnet[" << i << "] : ";
    for ( size_t j=0; j< cgr.bnets[i].candidate_routes.size(); j++ ) {
      double  max_decision_variable = 0;
      //  debug
      if ( debug_mode == true )
        cout << primalSol[primal_index+j] << " " << cgr.computeNetRouteWeight ( cgr.bnets[i].candidate_routes[j] ) << " ";
      //
      if ( primalSol[primal_index+j] > max_decision_variable ) {
        max_decision_variable = primalSol[primal_index+j];
        index_of_max_decision_variable = j;
      }
    }
    if ( debug_mode==true )
      cout << endl;
    assert ( index_of_max_decision_variable != -1 );
    primal_index += cgr.bnets[i].candidate_routes.size();
    //  update edges for bnets[i]
    cgr.bnets[i].edges.clear();
    cgr.bnets[i].edges = cgr.bnets[i].candidate_routes[index_of_max_decision_variable].edges;

    if ( debug_mode==true ) {
      if ( cgr.computeNetRouteWeight(cgr.bnets[i].candidate_routes[0]) != cgr.computeNetRouteWeight(cgr.bnets[i].candidate_routes[index_of_max_decision_variable])  ) {
        cout << "Bnet " << i<< "\t" << cgr.modified_net_type[i] << endl;
        cout << "\tcandidate rout " << 0 << "\t" << cgr.computeNetRouteWeight ( cgr.bnets[i].candidate_routes[0] ) << endl;
        cout << "\tcandidate rout " << 1 << "\t" << cgr.computeNetRouteWeight ( cgr.bnets[i].candidate_routes[1] ) << endl;
        cout << "\tcandidate rout " << index_of_max_decision_variable << "\t" << cgr.computeNetRouteWeight ( cgr.bnets[i].candidate_routes[index_of_max_decision_variable] ) << endl;
      }
    }
  }

  debug_mode = false;
  //  compute edge weight
  cgr.j_star.clear();
  for ( vector<int>::iterator itr_edge = cgr.E_id.begin(); itr_edge != cgr.E_id.end(); itr_edge++ ) {
    int best_j = 0;
    int current_dual_base_index = cgr.bnets.size() + (itr_edge - cgr.E_id.begin()) * cgr.jsize;
    //    if ( debug_mode==true ) 
    //      if ( cgr.E[*itr_edge].util>0 && cgr.E[*itr_edge].cap>0 && cgr.E[*itr_edge].congestion_ratio > cgr.critical_ratio )
    //        cgr.printEdge ( *itr_edge );
    for ( int k=0; k < cgr.jsize; k++ ) {
      if ( debug_mode == true )
        cout << dualSol[current_dual_base_index + k] << " ";
      if ( dualSol[current_dual_base_index + k] > dualSol[current_dual_base_index + best_j] )
        best_j = k;
    }
    cgr.j_star.push_back ( best_j );
    cgr.E[*itr_edge].pi = dualSol[current_dual_base_index + best_j ];
    cgr.E[*itr_edge].weight= cgr.mj[best_j] * cgr.E[*itr_edge].pi / (cgr.E[*itr_edge].cap+0.0001);
    if ( debug_mode == true )
      cout << " " << cgr.E[*itr_edge].weight << " | ";
  }
  delete [] primalSol;
  primalSol = NULL;
  delete [] reducedCost;
  reducedCost = NULL;

  return 0;
}

int Parser_cplex::CGR_AddNewRows (CPXENVptr env, CPXLPptr lp)
{
  //** allocate memory **

  //  cout << endl << "ADD LP ROWS" << endl;

  int status; 
  int rcnt = cgr.bnets.size() + cgr.E_id.size() * cgr.jsize;
  double * rhs = new double [rcnt];
  char * sense = new char [rcnt];

  int currIndex = 0;
  for (int i = 0; i < cgr.bnets.size(); i++)
  {
    rhs[currIndex] = 1.0;		
    sense[currIndex++] = 'E';			
  }

  for (int i = 0; i < cgr.E_id.size(); i++)
  {
    int edge_id = cgr.E_id[i];
    for ( int j=0; j < cgr.jsize; j++ )
    {
      rhs[currIndex] = cgr.bj[j] + cgr.mj[j] * cgr.E[edge_id].blk / ( cgr.E[edge_id].cap+0.0001);
      sense[currIndex++] = 'G';
    }
  }

  assert ( currIndex == rcnt );

  //** add rows to LP **

  status = CPXnewrows (env, lp, rcnt, rhs, sense, NULL, NULL);

  delete [] rhs;
  rhs = NULL;
  delete [] sense;
  sense = NULL;

  //  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  //  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int Parser_cplex::CGR_AddGridEdgeColumns (CPXENVptr env, CPXLPptr lp)
{
  //  **  VARS : p_e
  //      foreach e in E
  //
  //** allocate memory **
  int status = 1; 

  //  cout<<endl<<"ADD EDGE COLUMNS" << endl;

  for ( size_t i=0; i< cgr.E_id.size(); i++ ) {
    status = CGR_AddSingleEdge ( env, lp, i );
    assert ( status == 0 );
  }

  //  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  //  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int  Parser_cplex::CGR_AddSingleEdge ( CPXENVptr env, CPXLPptr lp, int edge ) 
{

  int edge_id = cgr.E_id[edge];

  bool  debug_mode = false;

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;

  int   ccnt = 1;
  int   nzcnt=  cgr.jsize;

  lb = new double[ccnt];    
  lb[0] = 0;
  obj = new double[ccnt];
  //  assert ( cgr.E[edge_id].type != "VIA" );
  //  assert ( cgr.E[edge_id].cap != 0 );
  obj[0] = 1;
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  for ( int i=0; i < cgr.jsize; i++ ) {
    cmatind[i] = cgr.bnets.size() + edge * cgr.jsize + i ;
    cmatval[i] = 1;
  }

  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg,
      cmatind, cmatval, lb, NULL, NULL );

  delete [] obj;
  obj = NULL;
  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;

  return  status;
}

int Parser_cplex::CGR_AddCoarseNetColumns (CPXENVptr env, CPXLPptr lp)
{
  //  cout<<endl<<"ADD NET COLUMNS" << endl;
  int status = 1;
  size_t   i,j;
  for ( i = 0; i < cgr.bnets.size(); i++ )  {
    for ( j=0; j < cgr.bnets[i].candidate_routes.size(); j++ ) {
      status = CGR_AddSingleRoute ( env, lp, make_pair(i,j) );
      assert ( status == 0 );
    }
  }

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;
  return (status);
}

int  Parser_cplex::CGR_AddSingleRoute ( CPXENVptr env, CPXLPptr lp, pair<int,int> xit ) 
{

  bool  debug_mode = false;

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;
  double * ub;

  int net_id = xit.first;
  int t = xit.second;
  NetRoute  &net_route = cgr.bnets[net_id].candidate_routes[t];
  /*  add a single column */
  int ccnt = 1;
  //  cnt route edges appeared in E_id 
  vector<int> edges;
  for ( set<int>::iterator it = net_route.edges.begin(); it!= net_route.edges.end(); it++ ) {
    if ( cgr.E[*it].cap==0 || cgr.E[*it].type=="VIA" )
      continue;
    edges.push_back (*it);
  }

  int lpRows = cgr.bnets.size() + edges.size() * cgr.jsize;

  /*  determine the number of non-zero entries */
  int nzcnt = 1 + edges.size() * cgr.jsize;
  lb = new double[ccnt];    
  lb[0] = 0;
  ub = new double[ccnt];
  ub[0] = 1;
  obj = new double[ccnt];
  obj[0] = 0;
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  cmatind[0] = net_id;
  cmatval[0] = 1;

  for ( size_t i=0; i< edges.size(); i++ ) {
    assert ( cgr.map_from_edge_id_to_index_in_E_id.count(edges[i]) == 1 );
    int index_in_E_id = cgr.map_from_edge_id_to_index_in_E_id[edges[i]];
    if ( debug_mode )  {
      cgr.printEdge ( edges[i] );
      cout << "after mapping : " << index_in_E_id << " in " << cgr.E_id.size() << endl;
    }
    for ( int j=0; j < cgr.jsize; j++ ) {
      cmatind[1 + i * cgr.jsize + j ] = cgr.bnets.size() + index_in_E_id * cgr.jsize + j;
      cmatval[1 + i * cgr.jsize + j ] = -cgr.mj[j] / ( cgr.E[edges[i]].cap + 0.0001 ) ;
    }
  }

  if ( debug_mode )
    cout << "net route (" << net_id << ","<< t << ").edges.size() = " << edges.size() << endl;

  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg,
      cmatind, cmatval, lb, ub, NULL );

  if ( status != 0 )

    if ( debug_mode ) {
      cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) <<endl;
      cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;
    }

  //** add new columns to LP **

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;

  return (status);
}

//---------------------------------------------------------------------------  
void Parser_cplex::Free_and_Null (char **ptr)
{
  if ( *ptr != NULL ) 
  {
    free (*ptr);
    *ptr = NULL;
  }
}
//---------------------------------------------------------------------------  

void Parser_cplex::Sleep (int sleeptime)
{
  int sec10 = sleeptime / 10;
  clock_t endsleep;

  for (int i = 0; i < sec10; i++)
  {
    endsleep = clock () + 10 * CLOCKS_PER_SEC ;
    while (clock() < endsleep)
    {

    }
  }
}
//---------------------------------------------------------------------------  

bool Parser_cplex::GetSelfLicense ()
{
  FILE * fp;
  int tmp;
  fp = fopen ("/pong/usr3/t/taihsuan/routing/license/cplex_license.txt", "r");

  if (fp == NULL)
  {
    fp = fopen ("/pong/usr3/t/taihsuan/routing/license/cplex_license.txt", "w");
    fprintf (fp, "1\n");
    fprintf (stderr, "cplex used license now : 1\n");
    fclose (fp);
    return true; 
  }	

  int licenseNow;
  tmp = fscanf (fp, "%d\n", &licenseNow);
  fclose (fp);

  if (licenseNow >= 7)
    return false;
  else
  {
    licenseNow += 1;

    fp = fopen ("/pong/usr3/t/taihsuan/routing/license/cplex_license.txt", "w");
    fprintf (fp, "%d\n", licenseNow);
    fprintf (stderr, "cplex used license now : %d\n", licenseNow);
    fclose (fp);
    return true; 		
  }
}
//---------------------------------------------------------------------------  

void Parser_cplex::ReleaseSelfLicense ()
{
  FILE * fp;
  int tmp;

  fp = fopen ("/pong/usr3/t/taihsuan/routing/license/cplex_license.txt", "r");

  if (fp == NULL)
  {
    return; 
  }	

  int licenseNow;
  tmp = fscanf (fp, "%d\n", &licenseNow);
  fclose (fp);

  licenseNow -= 1;

  if (licenseNow < 0 )
    licenseNow = 0;

  fp = fopen ("/pong/usr3/t/taihsuan/routing/license/cplex_license.txt", "w");
  fprintf (fp, "%d\n", licenseNow);
  fprintf (stderr, "cplex used license now : %d\n", licenseNow);
  fclose (fp);	
}

//  Added on NOV 6
int Parser_cplex::CreateLargeILP ()
{	
  cout<<endl<<"create large LP problem"<<endl << endl;

  int status = 0;

  CPXchgobjsen ( lpEnv, lp, CPX_MIN );

  cout << "# mini nets = " << cgr.mini_net_id.size() << endl;

  status = NewRows (lpEnv, lp);
  assert ( !status );

  status = AddEdgeColumns (lpEnv, lp);
  assert ( !status );

  status = AddCoarseNetColumns (lpEnv, lp);
  assert ( !status );

  status = SolveLP (lpEnv, lp);
  assert ( !status );

TERMINATE:	

  return (status);
}

int Parser_cplex::NewRows (CPXENVptr env, CPXLPptr lp)
{
  //** allocate memory **

  //  cout << endl << "ADD LP ROWS" << endl;

  int status; 
  int rcnt = cgr.mini_net_id.size() + cgr.E.size() * cgr.jsize;
  double * rhs = new double [rcnt];
  char * sense = new char [rcnt];

  int currIndex = 0;
  for (int i = 0; i < cgr.mini_net_id.size(); i++)
  {
    rhs[currIndex] = 1.0;		
    sense[currIndex++] = 'E';			
  }

  for (int i = 0; i < cgr.E.size(); i++)
  {
    for ( int j=0; j < cgr.jsize; j++ )
    {
      rhs[currIndex] = cgr.bj[j] + cgr.mj[j] * (cgr.E[i].blk+cgr.E[i].virtual_util) / ( cgr.E[i].cap+0.0001);
      sense[currIndex++] = 'G';
    }
  }

  assert ( currIndex == rcnt );

  //** add rows to LP **

  status = CPXnewrows (env, lp, rcnt, rhs, sense, NULL, NULL);

  delete [] rhs;
  rhs = NULL;
  delete [] sense;
  sense = NULL;

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int Parser_cplex::AddEdgeColumns (CPXENVptr env, CPXLPptr lp)
{
  //  **  VARS : p_e
  //      foreach e in E
  //
  //** allocate memory **
  int status = 1; 

  //  cout<<endl<<"ADD EDGE COLUMNS" << endl;

  for ( size_t i=0; i< cgr.E.size(); i++ ) {
    status = AddSingleEdge ( env, lp, i );
    assert ( status == 0 );
  }

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int  Parser_cplex::AddSingleEdge ( CPXENVptr env, CPXLPptr lp, int edge_id ) 
{
  bool  debug_mode = false;

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;

  int   ccnt = 1;
  int   nzcnt=  cgr.jsize;

  lb = new double[ccnt];    
  lb[0] = 0;
  obj = new double[ccnt];
  obj[0] = ( cgr.E[edge_id].cap!=0 );
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  for ( int i=0; i < cgr.jsize; i++ ) {
    cmatind[i] = cgr.mini_net_id.size() + edge_id* cgr.jsize + i ;
    cmatval[i] = 1;
  }

  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg,
      cmatind, cmatval, lb, NULL, NULL );

  delete [] obj;
  obj = NULL;
  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;

  return  status;
}

int Parser_cplex::AddCoarseNetColumns (CPXENVptr env, CPXLPptr lp)
{
  //  cout<<endl<<"ADD NET COLUMNS" << endl;
  int status = 1;
  size_t   i,j;
  for ( i = 0; i < cgr.mini_net_id.size(); i++ )  {
    assert ( cgr.nets[cgr.mini_net_id[i]].candidate_routes.size()>0 );
    for ( j=0; j < cgr.nets[cgr.mini_net_id[i]].candidate_routes.size(); j++ ) {
      status = AddSingleRoute ( env, lp, i, j );
      assert ( status == 0 );
    }
  }

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;
  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int  Parser_cplex::AddSingleRoute ( CPXENVptr env, CPXLPptr lp, int mini_id, int candidate_route_id ) 
{

  bool  debug_mode = false;

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;
  double * ub;

  /*  add a single column */
  int ccnt = 1;
  int lpRows = cgr.mini_net_id.size() + cgr.E.size() * cgr.jsize;
  /*  determine the number of non-zero entries */
  set<int> edges =cgr.nets[cgr.mini_net_id[mini_id]].candidate_routes[candidate_route_id].edges;
  int nzcnt = 1 + edges.size() * cgr.jsize;
  lb = new double[ccnt];    
  lb[0] = 0;
  ub = new double[ccnt];
  ub[0] = 1;
  obj = new double[ccnt];
  obj[0] = 0;
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  cmatind[0] = mini_id;
  cmatval[0] = 1;

  int curr_index = 1;
  for ( set<int>::iterator it = edges.begin(); it != edges.end(); it++ ) {
    for ( size_t j=0; j < cgr.jsize; j++ ) {
      cmatind[curr_index+j] = cgr.mini_net_id.size() + (*it)*cgr.jsize + j;
      cmatval[curr_index+j] = -cgr.mj[j] / ( cgr.E[*it].cap+0.0001 );
    }
    curr_index += cgr.jsize;
  }

  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg, cmatind, cmatval, lb, ub, NULL );

  //** add new columns to LP **

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;

  return (status);
}

int Parser_cplex::SolveLP (CPXENVptr env, CPXLPptr lp)
{
  int status;

  cout<<endl<<"solving LP problem"<<endl;

  /* compute the total number of candidate routes */
  cgr.x_var_size = 0;
  for ( vector<int>::iterator itr = cgr.mini_net_id.begin(); itr != cgr.mini_net_id.end(); itr ++ ) {
    cgr.x_var_size += cgr.nets[*itr].candidate_routes.size();
  }

  status = CPXsetdblparam (lpEnv, CPX_PARAM_TILIM, 10);
  if (status) {
    cerr << "Cannot set time limitation : " << status << endl;
    return  -1;
  }			

  status = CPXsetintparam (lpEnv, CPX_PARAM_LPMETHOD , CPX_ALG_BARRIER);
  if (status) {
    cerr << "Cannot set optimization method : " << status << endl;
    return  -1;
  }

  status = CPXlpopt (env, lp);
  if (status) {
    cerr << "Cannot optimize LP problem : " << status << endl;
    return  -1;
  }

  primalSol = new double [cgr.E.size() + cgr.x_var_size];

  dualSol = new double [cgr.mini_net_id.size() + cgr.E.size() * cgr.jsize];

  numcols = CPXgetnumcols ( lpEnv, lp );
  numrows = CPXgetnumrows ( lpEnv, lp );

  cout << "\tnumcols = " << numcols << endl;
  cout << "\tnumrows = " << numrows << endl;
  status = CPXsolution (env, lp, &solStat, &objValue, 
      primalSol, dualSol, NULL, NULL );
  exit(0);

  bool  debug_mode = false;
  double obj;

  /* update net route decision variables */
  int primal_index = cgr.E.size();

  for ( size_t i=0; i< cgr.mini_net_id.size(); i++ ) {
    int index_of_max_decision_variable = -1;
    for ( size_t j=0; j< cgr.nets[cgr.mini_net_id[i]].candidate_routes.size(); j++ ) {
      double  max_decision_variable = 0;
      if ( primalSol[primal_index+j] > max_decision_variable ) {
        max_decision_variable = primalSol[primal_index+j];
        index_of_max_decision_variable = j;
      }
    }
    assert ( index_of_max_decision_variable != -1 );
    primal_index += cgr.nets[cgr.mini_net_id[i]].candidate_routes.size();
    //  update edges for bnets[i]
    cgr.nets[cgr.mini_net_id[i]].edges.clear();
    cgr.nets[cgr.mini_net_id[i]].edges = cgr.nets[cgr.mini_net_id[i]].candidate_routes[index_of_max_decision_variable].edges;
  }

  //  compute edge weight
  cgr.j_star.clear();
  for ( size_t i=0; i < cgr.E.size(); i++ ) {
    int j_star = 0;
    int current_dual_base_index = cgr.mini_net_id.size() + i * cgr.jsize;
    for ( int k=1; k < cgr.jsize; k++ )
      if ( dualSol[current_dual_base_index + k] > dualSol[current_dual_base_index + j_star] )
        j_star = k;
    cgr.j_star.push_back ( j_star );
    cgr.E[i].pi = dualSol[current_dual_base_index + j_star];
    cgr.E[i].weight= cgr.mj[j_star] * cgr.E[i].pi / (cgr.E[i].cap+0.0001);
  }
  delete [] primalSol;
  primalSol = NULL;
  delete [] reducedCost;
  reducedCost = NULL;

  return 0;
}

//  Added on NOV 7
int Parser_cplex::CreateLargeILP2 ()
{	
  cout<<endl<<"create large LP problem"<<endl << endl;

  int status = 0;

  CPXchgobjsen ( lpEnv, lp, CPX_MIN );

  cout << "small nets = " << cgr.bnets.size() << endl;

  status = NewRows2 (lpEnv, lp);
  assert ( !status );

  status = AddEdgeColumns2 (lpEnv, lp);
  assert ( !status );

  status = AddCoarseNetColumns2 (lpEnv, lp);
  assert ( !status );

  status = SolveLP2 (lpEnv, lp);
  assert ( !status );

TERMINATE:	

  return (status);
}

int Parser_cplex::NewRows2 (CPXENVptr env, CPXLPptr lp)
{
  //** allocate memory **

  //  cout << endl << "ADD LP ROWS" << endl;

  int status; 
  int rcnt = cgr.bnets.size() + cgr.E_id.size() * cgr.jsize;
  double * rhs = new double [rcnt];
  char * sense = new char [rcnt];

  int currIndex = 0;
  for (int i = 0; i < cgr.bnets.size(); i++)
  {
    rhs[currIndex] = 1.0;		
    sense[currIndex++] = 'E';			
  }

  for (int i = 0; i < cgr.E_id.size(); i++)
  {
    int edge_id = cgr.E_id[i];
    for ( int j=0; j < cgr.jsize; j++ )
    {
      rhs[currIndex] = cgr.bj[j] + cgr.mj[j] * (cgr.E[edge_id].blk+cgr.E[edge_id].util) / ( cgr.E[edge_id].cap+0.0001); //  blk+util is fixed
      sense[currIndex++] = 'G';
    }
  }

  assert ( currIndex == rcnt );

  //** add rows to LP **

  status = CPXnewrows (env, lp, rcnt, rhs, sense, NULL, NULL);

  delete [] rhs;
  rhs = NULL;
  delete [] sense;
  sense = NULL;

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int Parser_cplex::AddEdgeColumns2 (CPXENVptr env, CPXLPptr lp)
{
  //  **  VARS : p_e
  //      foreach e in E
  //
  //** allocate memory **
  int status = 1; 

  //  cout<<endl<<"ADD EDGE COLUMNS" << endl;

  for ( size_t i=0; i< cgr.E_id.size(); i++ ) {
    status = AddSingleEdge2 ( env, lp, i );
    assert ( status == 0 );
  }

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int  Parser_cplex::AddSingleEdge2 ( CPXENVptr env, CPXLPptr lp, int edge_id ) 
{
  bool  debug_mode = false;

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;

  int   ccnt = 1;
  int   nzcnt=  cgr.jsize;

  lb = new double[ccnt];    
  lb[0] = 0;
  obj = new double[ccnt];
  obj[0] = 1;
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  for ( int j=0; j < cgr.jsize; j++ ) {
    cmatind[j] = cgr.bnets.size() + edge_id* cgr.jsize + j ;
    cmatval[j] = 1;
  }

  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg,
      cmatind, cmatval, lb, NULL, NULL );

  delete [] obj;
  obj = NULL;
  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;

  return  status;
}

int Parser_cplex::AddCoarseNetColumns2 (CPXENVptr env, CPXLPptr lp)
{
  //  cout<<endl<<"ADD NET COLUMNS" << endl;
  int status = 1;
  size_t   i,j;
  for ( i = 0; i < cgr.bnets.size(); i++ )  {
    assert ( cgr.bnets[i].candidate_routes.size()>0 );
    for ( j=0; j < cgr.bnets[i].candidate_routes.size(); j++ ) {
      status = AddSingleRoute2 ( env, lp, i, j );
      assert ( status == 0 );
    }
  }

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;
  return (status);

}

int  Parser_cplex::AddSingleRoute2 ( CPXENVptr env, CPXLPptr lp, int bnet_id, int candidate_route_id ) 
{

  bool  debug_mode = false;

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;
  double * ub;

  int ccnt = 1;
  int lpRows = cgr.bnets.size() + cgr.E_id.size() * cgr.jsize;
  /*  determine the number of non-zero entries */
  NetRoute  &net_route = cgr.bnets[bnet_id].candidate_routes[candidate_route_id];
  int nzcnt = 1 + net_route.edges.size() * cgr.jsize;
  lb = new double[ccnt];    
  lb[0] = 0;
  ub = new double[ccnt];
  ub[0] = 1;
  obj = new double[ccnt];
  obj[0] = 0;
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  cmatind[0] = bnet_id;
  cmatval[0] = 1;

  vector<int> edges ( net_route.edges.begin(), net_route.edges.end() );
  for ( size_t i=0; i< edges.size(); i++ ) {
    assert ( cgr.map_from_edge_id_to_index_in_E_id.count(edges[i]) == 1 );
    int index_in_E_id = cgr.map_from_edge_id_to_index_in_E_id[edges[i]];
    for ( int j=0; j < cgr.jsize; j++ ) {
      cmatind[1 + i * cgr.jsize + j ] = cgr.bnets.size() + index_in_E_id * cgr.jsize + j;
      cmatval[1 + i * cgr.jsize + j ] = -cgr.mj[j] / ( cgr.E[edges[i] ].cap + 0.0001 ) ;
    }
  }
  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg,
      cmatind, cmatval, lb, ub, NULL );

  //** add new columns to LP **

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;

  return (status);
}

int Parser_cplex::SolveLP2 (CPXENVptr env, CPXLPptr lp)
{
  int status;

  cout<<endl<<"solving LP problem"<<endl;

  /* compute the total number of candidate routes */
  cgr.x_var_size = 0;
  for ( vector<CoarseNet>::iterator itr_net = cgr.bnets.begin(); itr_net != cgr.bnets.end(); itr_net ++ ) {
    cgr.x_var_size += (*itr_net).candidate_routes.size();
  }

  status = CPXsetdblparam (lpEnv, CPX_PARAM_TILIM, 10);
  if (status) {
    cerr << "Cannot set time limitation : " << status << endl;
    return  -1;
  }			

  status = CPXsetintparam (lpEnv, CPX_PARAM_LPMETHOD , CPX_ALG_BARRIER);
  if (status) {
    cerr << "Cannot set optimization method : " << status << endl;
    return  -1;
  }

  status = CPXlpopt (env, lp);
  if (status) {
    cerr << "Cannot optimize LP problem : " << status << endl;
    return  -1;
  }

  primalSol = new double [cgr.E_id.size() + cgr.x_var_size];

  dualSol = new double [cgr.bnets.size() + cgr.E_id.size() * cgr.jsize];

  numcols = CPXgetnumcols ( lpEnv, lp );
  numrows = CPXgetnumrows ( lpEnv, lp );

  status = CPXsolution (env, lp, &solStat, &objValue, 
      primalSol, dualSol, NULL, NULL );

  cout << "In this iteratoin : " << endl;
  cout << "\tnumcols = " << numcols << endl;
  cout << "\tnumrows = " << numrows << endl;
  cout << "\tobjval = " << objValue << endl;

  bool  debug_mode = false;
  double obj;

  /* update net route decision variables */
  int primal_index = cgr.E_id.size();

  for ( size_t i=0; i< cgr.bnets.size(); i++ ) {
    int index_of_max_decision_variable = -1;
    for ( size_t j=0; j< cgr.bnets[i].candidate_routes.size(); j++ ) {
      double  max_decision_variable = 0;
      if ( primalSol[primal_index+j] > max_decision_variable ) {
        max_decision_variable = primalSol[primal_index+j];
        index_of_max_decision_variable = j;
      }
    }
    assert ( index_of_max_decision_variable != -1 );
    primal_index += cgr.bnets[i].candidate_routes.size();
    //  update edges for bnets[i]
    cgr.bnets[i].edges.clear();
    cgr.bnets[i].edges = cgr.bnets[i].candidate_routes[index_of_max_decision_variable].edges;
  }

  //  compute edge weight
  cgr.j_star.clear();
  for ( vector<int>::iterator itr_edge = cgr.E_id.begin(); itr_edge != cgr.E_id.end(); itr_edge++ ) {
    int current_dual_base_index = cgr.bnets.size() + (itr_edge - cgr.E_id.begin()) * cgr.jsize;
    //  start computing edge weight
    cgr.E[*itr_edge].weight = 0;
    for ( int j=0; j < cgr.jsize; j++ )
      cgr.E[*itr_edge].weight += cgr.mj[j] * dualSol[current_dual_base_index + j] / ( cgr.E[*itr_edge].cap+0.0001 );
  }
  delete [] primalSol;
  primalSol = NULL;
  delete [] reducedCost;
  reducedCost = NULL;

  return 0;
}

//  Added on NOV 19 
int Parser_cplex::CreateLargeILP3 ( std::ofstream &stream_file )
{	
  cout<<endl<<"create large LP problem"<<endl << endl;

  int status = 0;

  CPXchgobjsen ( lpEnv, lp, CPX_MIN );

  status = NewRows3 (lpEnv, lp);
  assert ( !status );

  status = AddEdgeColumns3 (lpEnv, lp);
  assert ( !status );

  status = AddCoarseNetColumns3 (lpEnv, lp);
  assert ( !status );

  status = SolveLP3 (lpEnv, lp, stream_file );
  assert ( !status );

TERMINATE:	

  return (status);
}

int Parser_cplex::NewRows3 (CPXENVptr env, CPXLPptr lp)
{
  //** allocate memory **

  //  cout << endl << "ADD LP ROWS" << endl;

  int status; 
  int rcnt = cgr.bnets.size() + cgr.E_id.size() * cgr.jsize;
  double * rhs = new double [rcnt];
  char * sense = new char [rcnt];

  int currIndex = 0;
  for (int i = 0; i < cgr.bnets.size(); i++)
  {
    rhs[currIndex] = 1.0;		
    sense[currIndex++] = 'E';			
  }

  for (int i = 0; i < cgr.E_id.size(); i++)
  {
    int edge_id = cgr.E_id[i];
    for ( int j=0; j < cgr.jsize; j++ )
    {
      rhs[currIndex] = cgr.bj[j] + cgr.mj[j] * (cgr.E[edge_id].blk+cgr.E[edge_id].util) / ( cgr.E[edge_id].cap+0.0001); 
      sense[currIndex++] = 'G';
    }
  }

  assert ( currIndex == rcnt );

  //** add rows to LP **

  status = CPXnewrows (env, lp, rcnt, rhs, sense, NULL, NULL);

  delete [] rhs;
  rhs = NULL;
  delete [] sense;
  sense = NULL;

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int Parser_cplex::AddEdgeColumns3 (CPXENVptr env, CPXLPptr lp)
{
  //  **  VARS : p_e
  //      foreach e in E
  //
  //** allocate memory **
  int status = 1; 

  //  cout<<endl<<"ADD EDGE COLUMNS" << endl;

  for ( size_t i=0; i< cgr.E_id.size(); i++ ) {
    status = AddSingleEdge3 ( env, lp, i );
    assert ( status == 0 );
  }

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int  Parser_cplex::AddSingleEdge3 ( CPXENVptr env, CPXLPptr lp, int i ) 
{
  bool  debug_mode = false;
  int edge_id = cgr.E_id[i];

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;

  int   ccnt = 1;
  int   nzcnt=  cgr.jsize;

  lb = new double[ccnt];    
  lb[0] = 0;
  obj = new double[ccnt];
  obj[0] = (cgr.E[edge_id].cap!=0&&cgr.E[edge_id].type!="VIA")?1:0;
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  for ( int j=0; j < cgr.jsize; j++ ) {
    cmatind[j] = cgr.bnets.size() + i * cgr.jsize + j ;
    cmatval[j] = 1;
  }

  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg,
      cmatind, cmatval, lb, NULL, NULL );

  delete [] obj;
  obj = NULL;
  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;

  return  status;
}

int Parser_cplex::AddCoarseNetColumns3 (CPXENVptr env, CPXLPptr lp)
{
  //  cout<<endl<<"ADD NET COLUMNS" << endl;
  int status = 1;
  size_t   i,j;
  for ( i = 0; i < cgr.bnets.size(); i++ )  {
    assert ( cgr.bnets[i].candidate_routes.size()>0 );
    for ( j=0; j < cgr.bnets[i].candidate_routes.size(); j++ ) {
      status = AddSingleRoute3 ( env, lp, i, j );
      assert ( status == 0 );
    }
  }

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;
  return (status);

}

int  Parser_cplex::AddSingleRoute3 ( CPXENVptr env, CPXLPptr lp, int bnet_id, int candidate_route_id ) 
{

  bool  debug_mode = false;

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;
  double * ub;

  int ccnt = 1;
  /*  determine the number of non-zero entries */
  NetRoute  &net_route = cgr.bnets[bnet_id].candidate_routes[candidate_route_id];
  int nzcnt = 1 + net_route.edges.size() * cgr.jsize;
  lb = new double[ccnt];    
  lb[0] = 0;
  ub = new double[ccnt];
  ub[0] = 1;
  obj = new double[ccnt];
  obj[0] = 0;
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  cmatind[0] = bnet_id;
  cmatval[0] = 1;

  vector<int> edges ( net_route.edges.begin(), net_route.edges.end() );
  for ( size_t i=0; i< edges.size(); i++ ) {
    assert ( cgr.map_from_edge_id_to_index_in_E_id.count(edges[i]) == 1 );
    int index_in_E_id = cgr.map_from_edge_id_to_index_in_E_id[edges[i]];
    for ( int j=0; j < cgr.jsize; j++ ) {
      cmatind[1 + i * cgr.jsize + j ] = cgr.bnets.size() + index_in_E_id * cgr.jsize + j;
      cmatval[1 + i * cgr.jsize + j ] = -cgr.mj[j] / ( cgr.E[edges[i] ].cap + 0.0001 ) ;
    }
  }
  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg,
      cmatind, cmatval, lb, ub, NULL );

  //** add new columns to LP **

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;

  return (status);
}

int Parser_cplex::SolveLP3 (CPXENVptr env, CPXLPptr lp, std::ofstream &stream_file )
{
  int status;

  cout<<endl<<"solving LP problem"<<endl;

  /* compute the total number of candidate routes */
  cgr.x_var_size = 0;
  for ( vector<CoarseNet>::iterator itr_net = cgr.bnets.begin(); itr_net != cgr.bnets.end(); itr_net ++ ) {
    cgr.x_var_size += (*itr_net).candidate_routes.size();
  }

  status = CPXsetdblparam (lpEnv, CPX_PARAM_TILIM, 10);
  if (status) {
    cerr << "Cannot set time limitation : " << status << endl;
    return  -1;
  }			

  status = CPXsetintparam (lpEnv, CPX_PARAM_LPMETHOD , CPX_ALG_BARRIER);
  if (status) {
    cerr << "Cannot set optimization method : " << status << endl;
    return  -1;
  }

  status = CPXlpopt (env, lp);
  if (status) {
    cerr << "Cannot optimize LP problem : " << status << endl;
    return  -1;
  }

  primalSol = new double [cgr.E_id.size() + cgr.x_var_size];
  cout << "x_var_size() = " << cgr.x_var_size << endl;

  dualSol = new double [cgr.bnets.size() + cgr.E_id.size() * cgr.jsize];
  cout << "bnets.size()=" << cgr.bnets.size() << " E_id.size()=" << cgr.E_id.size() << " jsize=" << cgr.jsize << endl;

  numcols = CPXgetnumcols ( lpEnv, lp );
  numrows = CPXgetnumrows ( lpEnv, lp );

  status = CPXsolution (env, lp, &solStat, &objValue, 
      primalSol, dualSol, NULL, NULL );

  cout << "In this iteration: " << endl;
  cout << "\tnumcols = " << numcols << endl;
  cout << "\tnumrows = " << numrows << endl;
  cout << "\tobjval = " << objValue << endl;

  bool  debug_mode = false;
  double obj;

  /* update net route decision variables */
  int primal_index = cgr.E_id.size();

  cgr.consider_continuous = true;

  if ( cgr.consider_continuous == true ) {
    //  compute the continuous thing
    for ( vector<int>::iterator itr = cgr.E_id.begin(); itr!= cgr.E_id.end(); itr++ )
      cgr.E[*itr].virtual_util_fp = (double) cgr.E[*itr].util;
    int curr_index = cgr.E_id.size();
    for ( size_t i=0; i<cgr.bnets.size(); i++ ) {
      for ( size_t j=0; j< cgr.bnets[i].candidate_routes.size(); j++ ) {
        //  xit w\ curr_index
        NetRoute &route = cgr.bnets[i].candidate_routes[j];
        for ( set<int>::iterator it = route.edges.begin(); it!= route.edges.end(); it++ ) {
          cgr.E[*it].virtual_util_fp  = cgr.E[*it].virtual_util_fp + primalSol[curr_index];
        }
        ++ curr_index;
      }
    }
    cout << "2" << endl;
    assert ( curr_index == numcols );

    double  p1 = 0;
    double  p2 = 0;
    double  p3 = 0;
    double  p4 = 0;
    int n1 = 0;
    int n2 = 0;
    int n3 = 0;
    int n4 = 0;

    ofstream cr1;
    cr1.open( "cr1.txt" );

    cr1 << "helloworld" << endl;

    int tot_cnt =0;

    for ( vector<int>::iterator itr = cgr.E_id.begin(); itr!= cgr.E_id.end(); itr++ ) {
      double  tmp_cr = (double)( cgr.E[*itr].blk+cgr.E[*itr].virtual_util_fp ) / ( cgr.E[*itr].cap + 0.0001 );
      //  compute distribution
      if ( cgr.E[*itr].cap == 0 || cgr.E[*itr].type == "VIA" )
        continue;
      cr1 << tmp_cr << endl;
      ++ tot_cnt;
      if ( tmp_cr < cgr.bdr1) {
        p1 += cgr.getPenalty ( tmp_cr );
        n1++;
      }
      else  if ( tmp_cr < cgr.bdr2) {
        p2 += cgr.getPenalty ( tmp_cr );
        n2++;
      }
      else  if ( tmp_cr < cgr.bdr3 ) {
        p3 += cgr.getPenalty ( tmp_cr );
        n3++;
      }
      else {
        p4 += cgr.getPenalty ( tmp_cr );
        n4++;
      }
      cgr.E[*itr].congestion_ratio  = tmp_cr;   //  update the edge congestion ratio using continuous thing
    }
    cr1 << "EOF" << endl;

    cout << "## Continuous distribution" << endl;
    cout << n1 << " " << n2 << " " << n3 << " " << n4 << " ";
    cout << p1/(p1+p2+p3+p4) << " "; 
    cout << p2/(p1+p2+p3+p4) << " "; 
    cout << p3/(p1+p2+p3+p4) << " "; 
    cout << p4/(p1+p2+p3+p4) << endl;

    stream_file << "## Continuous distribution" << endl;
    stream_file << n1 << " " << n2 << " " << n3 << " " << n4 << endl;
    stream_file << p1/(p1+p2+p3+p4) << " "; 
    stream_file << p2/(p1+p2+p3+p4) << " "; 
    stream_file << p3/(p1+p2+p3+p4) << " "; 
    stream_file << p4/(p1+p2+p3+p4) << endl;
  }

  cout << "update net info" << endl;
  for ( size_t i=0; i< cgr.bnets.size(); i++ ) {
    int index_of_max_decision_variable = -1;
    for ( size_t j=0; j< cgr.bnets[i].candidate_routes.size(); j++ ) {
      double  max_decision_variable = 0;
      if ( primalSol[primal_index+j] > max_decision_variable ) {
        max_decision_variable = primalSol[primal_index+j];
        index_of_max_decision_variable = j;
      }
    }
    assert ( index_of_max_decision_variable != -1 );
    primal_index += cgr.bnets[i].candidate_routes.size();
    //  update edges for bnets[i]
    cgr.bnets[i].edges.clear();
    cgr.bnets[i].edges = cgr.bnets[i].candidate_routes[index_of_max_decision_variable].edges;
  }

  cout << "update edge info" << endl;
  //
  //  compute edge weight
  cgr.j_star.clear();
  for ( vector<int>::iterator itr_edge = cgr.E_id.begin(); itr_edge != cgr.E_id.end(); itr_edge++ ) {
    int current_dual_base_index = cgr.bnets.size() + (itr_edge - cgr.E_id.begin()) * cgr.jsize;
    //  start computing edge weight
    cgr.E[*itr_edge].weight = 0;
    for ( int j=0; j < cgr.jsize; j++ )
      cgr.E[*itr_edge].weight += cgr.mj[j] * dualSol[current_dual_base_index + j] / ( cgr.E[*itr_edge].cap+0.0001 );
  }

  cout << "done" << endl;
  delete [] primalSol;
  primalSol = NULL;
  delete [] reducedCost;
  reducedCost = NULL;

  return 0;
}

//  Added on NOV 27 
int Parser_cplex::CreateLargeILP4 ()
{	
  cout<<endl<<"create large LP problem"<<endl << endl;

  int status = 0;

  CPXchgobjsen ( lpEnv, lp, CPX_MIN );

  status = NewRows4 (lpEnv, lp);
  assert ( !status );

  status = AddEdgeColumns4 (lpEnv, lp);
  assert ( !status );

  status = AddCoarseNetColumns4 (lpEnv, lp);
  assert ( !status );

  status = SolveLP4 (lpEnv, lp);
  assert ( !status );

TERMINATE:	

  return (status);
}

int Parser_cplex::NewRows4 (CPXENVptr env, CPXLPptr lp)
{
  //** allocate memory **

  //  cout << endl << "ADD LP ROWS" << endl;

  int status; 
  int rcnt = cgr.bnets.size() + cgr.E_id.size();
  double * rhs = new double [rcnt];
  char * sense = new char [rcnt];

  int currIndex = 0;
  for (int i = 0; i < cgr.bnets.size(); i++)
  {
    rhs[currIndex] = 1.0;		
    sense[currIndex++] = 'E';			
  }

  for (int i = 0; i < cgr.E_id.size(); i++)
  {
    int edge_id = cgr.E_id[i];
    rhs[currIndex] = cgr.slope * ((cgr.E[edge_id].blk+cgr.E[edge_id].util) / ( cgr.E[edge_id].cap+0.0001) - cgr.TCR); 
    //  blk+util is fixed
    sense[currIndex++] = 'G';
  }

  assert ( currIndex == rcnt );

  //** add rows to LP **

  status = CPXnewrows (env, lp, rcnt, rhs, sense, NULL, NULL);

  delete [] rhs;
  rhs = NULL;
  delete [] sense;
  sense = NULL;

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int Parser_cplex::AddEdgeColumns4 (CPXENVptr env, CPXLPptr lp)
{
  //  **  VARS : p_e
  //      foreach e in E
  //
  //** allocate memory **
  int status = 1; 

  //  cout<<endl<<"ADD EDGE COLUMNS" << endl;

  for ( size_t i=0; i< cgr.E_id.size(); i++ ) {
    status = AddSingleEdge4 ( env, lp, i );
    assert ( status == 0 );
  }

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  return (status);
}

int  Parser_cplex::AddSingleEdge4 ( CPXENVptr env, CPXLPptr lp, int i ) 
{
  bool  debug_mode = false;
  int edge_id = cgr.E_id[i];

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;

  int   ccnt = 1;
  int   nzcnt=  1;

  lb = new double[ccnt];    
  lb[0] = 0;
  obj = new double[ccnt];
  obj[0] = (cgr.E[edge_id].cap!=0&&cgr.E[edge_id].type!="VIA")?1:0;
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  cmatind[0] = cgr.bnets.size() + i ;
  cmatval[0] = 1;

  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg,
      cmatind, cmatval, lb, NULL, NULL );

  delete [] obj;
  obj = NULL;
  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;

  return  status;
}

int Parser_cplex::AddCoarseNetColumns4 (CPXENVptr env, CPXLPptr lp)
{
  //  cout<<endl<<"ADD NET COLUMNS" << endl;
  int status = 1;
  size_t   i,j;
  for ( i = 0; i < cgr.bnets.size(); i++ )  {
    assert ( cgr.bnets[i].candidate_routes.size()>0 );
    for ( j=0; j < cgr.bnets[i].candidate_routes.size(); j++ ) {
      status = AddSingleRoute4 ( env, lp, i, j );
      assert ( status == 0 );
    }
  }

  cout << "numcols = " << CPXgetnumcols ( lpEnv, lp ) << endl;
  cout << "numrows = " << CPXgetnumrows ( lpEnv, lp ) << endl;

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;
  return (status);

}

int  Parser_cplex::AddSingleRoute4 ( CPXENVptr env, CPXLPptr lp, int bnet_id, int candidate_route_id ) 
{

  bool  debug_mode = false;

  int  status;
  double * obj;
  int * cmatbeg;
  int * cmatind;
  double * cmatval;
  double * lb;
  double * ub;

  int ccnt = 1;
  /*  determine the number of non-zero entries */
  NetRoute  &net_route = cgr.bnets[bnet_id].candidate_routes[candidate_route_id];
  int nzcnt = 1 + net_route.edges.size();
  lb = new double[ccnt];    
  lb[0] = 0;
  ub = new double[ccnt];
  ub[0] = 1;
  obj = new double[ccnt];
  obj[0] = 0;
  cmatbeg = new int [ccnt];
  cmatbeg[0] = 0;

  cmatind = new int [nzcnt];
  cmatval = new double [nzcnt];

  cmatind[0] = bnet_id;
  cmatval[0] = 1;

  vector<int> edges ( net_route.edges.begin(), net_route.edges.end() );
  for ( size_t i=0; i< edges.size(); i++ ) {
    assert ( cgr.map_from_edge_id_to_index_in_E_id.count(edges[i]) == 1 );
    int index_in_E_id = cgr.map_from_edge_id_to_index_in_E_id[edges[i]];
    cmatind[1 + i ] = cgr.bnets.size() + index_in_E_id;
    cmatval[1 + i ] = - cgr.slope / ( cgr.E[edges[i] ].cap + 0.0001 ) ;
  }
  status = CPXaddcols ( env, lp, ccnt, nzcnt, obj, cmatbeg,
      cmatind, cmatval, lb, ub, NULL );

  //** add new columns to LP **

  delete [] cmatbeg;
  cmatbeg = NULL;
  delete [] cmatind;
  cmatind = NULL;
  delete [] cmatval;
  cmatval = NULL;
  delete [] lb;
  lb = NULL;
  delete [] ub;
  ub = NULL;

  return (status);
}

int Parser_cplex::SolveLP4 (CPXENVptr env, CPXLPptr lp)
{
  int status;

  cout<<endl<<"solving LP problem"<<endl;

  /* compute the total number of candidate routes */
  cgr.x_var_size = 0;
  for ( vector<CoarseNet>::iterator itr_net = cgr.bnets.begin(); itr_net != cgr.bnets.end(); itr_net ++ ) {
    cgr.x_var_size += (*itr_net).candidate_routes.size();
  }

  status = CPXsetdblparam (lpEnv, CPX_PARAM_TILIM, 10);
  if (status) {
    cerr << "Cannot set time limitation : " << status << endl;
    return  -1;
  }			

  status = CPXsetintparam (lpEnv, CPX_PARAM_LPMETHOD , CPX_ALG_BARRIER);
  if (status) {
    cerr << "Cannot set optimization method : " << status << endl;
    return  -1;
  }

  status = CPXlpopt (env, lp);
  if (status) {
    cerr << "Cannot optimize LP problem : " << status << endl;
    return  -1;
  }

  primalSol = new double [cgr.E_id.size() + cgr.x_var_size];
  cout << "x_var_size() = " << cgr.x_var_size << endl;

  dualSol = new double [cgr.bnets.size() + cgr.E_id.size()];
  cout << "bnets.size()=" << cgr.bnets.size() << " E_id.size()=" << cgr.E_id.size() << endl;

  numcols = CPXgetnumcols ( lpEnv, lp );
  numrows = CPXgetnumrows ( lpEnv, lp );

  status = CPXsolution (env, lp, &solStat, &objValue, 
      primalSol, dualSol, NULL, NULL );

  cout << "In this iteration: " << endl;
  cout << "\tnumcols = " << numcols << endl;
  cout << "\tnumrows = " << numrows << endl;
  cout << "\tobjval = " << objValue << endl;

  bool  debug_mode = false;
  double obj;

  /* update net route decision variables */
  int primal_index = cgr.E_id.size();

  cout << "update net info" << endl;
  for ( size_t i=0; i< cgr.bnets.size(); i++ ) {
    int index_of_max_decision_variable = -1;
    for ( size_t j=0; j< cgr.bnets[i].candidate_routes.size(); j++ ) {
      double  max_decision_variable = 0;
      if ( primalSol[primal_index+j] > max_decision_variable ) {
        //  debug
//        cout << primalSol[primal_index+j] << " ";
        max_decision_variable = primalSol[primal_index+j];
        index_of_max_decision_variable = j;
      }
//      cout << endl;
    }
    assert ( index_of_max_decision_variable != -1 );
    primal_index += cgr.bnets[i].candidate_routes.size();
    //  update edges for bnets[i]
    cgr.bnets[i].edges.clear();
    cgr.bnets[i].edges = cgr.bnets[i].candidate_routes[index_of_max_decision_variable].edges;
  }

  cout << "update edge info" << endl;
  //
  //  compute edge weight
  cgr.j_star.clear();
  for ( vector<int>::iterator itr_edge = cgr.E_id.begin(); itr_edge != cgr.E_id.end(); itr_edge++ ) 
    cgr.E[*itr_edge].weight = cgr.slope * dualSol[cgr.bnets.size()+itr_edge-cgr.E_id.begin()] / ( cgr.E[*itr_edge].cap+0.0001 );

  cout << "done" << endl;
  delete [] primalSol;
  primalSol = NULL;
  delete [] reducedCost;
  reducedCost = NULL;

  return 0;
}
